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Abstract. Here's a LISP library ot mathematical functions that calculate hyperbolic and inverse 
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functions. Mukiple linear regression. Fletcher-Powell unconstrained minimization, numerical 
integration, root finding, and convergence. Code to factor numbers and to do the Solovay-Strassen 
probabilistic prime test. 
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1. Introduction 

Several times in science and engineering we need to evaluate rather exotic functions -- elliptic 
integrals, Bessel functions, and probability density flmctions, for example. Here's where you can get 
your hands on some LISP code tliat will evaluate some of tliese functions. The code runs in both 
MACLISP and LISPM LISP. 

The functions in the library are scattered in several different files. The descriptions below will 
describe tlie function name, how it should be used, and the file where that function can be obtained. For 
example, die function FACTOR is in PS:<GLR.FUNCT>MOD.LSP. Before it can be used, Uie code 
must first be loaded by doing one of: 

(LOAD "PS:<GLR.FUNCT>MOD.FASL") ; For MACLISP 

(LOAD "OZ:PS:<GLR.FUNCT>MOD.QFASL") ; For MIT LISPM 

(LOAD "OZ:PS:<GLR.FUNCT>MOD.QBIN") ; For Symbolics LM2 

(LOAD "OZ:PS:<GLR.FUNCT>MOD.BIN") ; For Symbolics 3600 

After tlie file has been loaded, just call die function with the right arguments: 

(FACTOR 12345678987654321) ==> (3 3 3 3 37 37 333667 333667) 

If you have to compile one of the library files or you want to read the source version of a library 
file into your LISP, tlien you will have to load a special macro called IMPORT-FILE that is defined in 
tlie file PS:<GLR.LISP>IMPORT.LSP. You may want to use tlie IMPORT-FILE macro yourself. Its 
fomiat is: 

(IMPORT-FILE "OZ:PS:<GLR.FUNCT>MOD.LSP") 

The argument is die source filename wiUi a LISPM style hostname (the MACLISP version of tlie macro 
knows to ignore the host). The behavior of tlie macro depends on whether it is being compiled, 
evaluated, or loaded. If the IMPORT-FILE fomi is evaluated or loaded (EVAL-WHEN (EVAL LOAD) 
...) and tlie imported file has not already been loaded, it will load a binary file if it can find one or the 
source file if it cannot. At compile time (EVAL-WHEN (COMPILE) ...), die IMPORT-FILE forni turns 
into 

(PROGN 'COMPILE <first form in source fi1e> <load code>) 

By convention, die first form in die imported source file should be a DECLARE of all die external 
finictions and variables that die imported file defines. Including die declarations allows some compile 
time error checking and eliminates spurious warnings about undefined fimctions. llie <Ioadcode> 
conspires to load die binary or source file when when die file we are currently compiling gets loaded. 



2 NUMBKR TlinORY AND COMBINATORIAL FUNCTIONS 



2. Number Theory and Combinatorial Functions 

(FACTORIAL N) ps:<glr.funct>combin.lsp 

The number of ways to permute N objects. 

(CHOOSE N M) PS:<GLR.FUNCT>COMBIN.LSP 

The number of ways to place M distinguishable objects on N distinguishable plates. 

(BELL-NUMBER N) ps:<glr.funct>combin.lsp 

The number of ways to place N distinguishable objects on N indistinguishable plates. 

(CATALAN-NUMBER N) ps:<glr.funct>combin.lsp 

The number of ways to fully parendiesize a string of n symbols. 

(FIB N) PS:<GLR.FUNCT>FIB.LSP 

The Nth (FIXNUM) Fibonacci number. FIB(0) = 0, FIB(1) = 1, .... Uses a log(N) algoriOim. 

(FACTOR N) PS:<GLR.FUNCT>MOD.LSP 

Find the prime factors of die integer N (FIXNUM or BIGNUM). Returns a sorted list of die 
factors. 

(PRIME-TEST N &optional (TRIALS 50.)) ps:<glr.funct>mod.lsp 

Tests die primality of N using a probabilistic algoridim (see <Solovay>). Returns NIL, PRIME, 
or PROBABLY-PRIME (P{error} = 2-'ri^JALS) trials j^^st be FIXNUM.. 

(SMALLER-PRIME N) ps:<glr.funct>mod.lsp 

Finds die first prime diat is smaller dian N. Repeatedly calls PRIME-TEST. 

( JACOBI-SYMBOL P Q) ps:<glr.funct>mod.lsp 

Computes the Jacobi Symbol of P and Q. 

(TOTIENT N) ps:<glr.funct>mod.lsp 

Computes Euler's toUent function of N. Calls FACTOR as a subroudne. 



3. Discrete Fourier Transform 

These functions are used to compute discrete fourier transforms <Oppenheimer>. The inputs 
are two FLONUM arrays that represent die real and imaginary parts of die input. I'he algoridims require 
input arrays for die transform to be bit reversed; there are functions for doing die reversal, llic functions 
also require some inidalizadon. The length of die DPT must be a power of 2. 



3 DISCRiniv I OURIF.R FRANSFORM 



XIk] = 2 ^'2o exp(-j2 7rkm/N) 
x[k] = (1/N) 2 ^^'2q cxp( j 2 TT k m / N) 

(DFT-INIT LOGN) ps:<glr.funct>dft.lsp 

DFF-INrr returns a structure containing some tables that are needed by the OFF algorithm. 
This stnicture should be saved away and given to the DFf routine each time it is called. The structure 
need only be computed once for each size DFT. The length of tlie DFF is 2 ^^^ . 

(DFT-FORWARD X-REAL X-IMAG TABLES) ps:<glr.funct>dft.lsp 

This function does a discrete Fourier ti"ansform. The results are written back into X-REAL and 
XTMAG (the previous contents are lost) and are in sequential order (ie -- not bit reversed). TABLES is 
as returned by DPTTNIT. 

(DFT-REVERSE X-REAL X-IMAG TABLES) PS:<GLR.FUNCT>DFr.LSP 

This is just like DPT- FOR WARD, but it does the inverse transform. The input X arrays should 
be in sequential (not bit-reversed) order. The resulting arrays are also in sequential order. 

(DFT-REVERSE-ARRAY ARRAY TABLES) ps:<glr.funct>dft.lsp 

This bit reverses the first 2^^^^^ elements of ARRAY. 

(DFT-610 X-REAL X-IMAG TABLES) ps:<glr.funct>dft.lsp 

This function does a forward transform. X-REAL and X-IMAG are bit-reversed input arrays 
(to reverse them, use DFT-REVERSE-ARRAY). The results are returned in X-REAL and X-IMAG in 
sequential order. 

(DFT-618 X-REAL X-IMAG TABLES) ps:<glr.funct>dft.lsp 

The reverse transfonn. Input arrays are sequential order, output arrays arc bit-reversed. Use 
DFT-REVERSE-ARRAY to put them in sequential order. 



4. Trigonometric and other Functions 

These approximations come from <Abramowitz>, <DEC>, and <Hastings>. 

(EXPIO X) PS:<GLR.FUNCT>EXTFCN.LSP 

Computes lOtX. 



4 IRIGONOMEI RIC AND OJHliR I-UNCTIONS 



(LOGIO X) 

Log base lOofX. 

(TAN X) 
(SEC X) 
(CSC X) 
(COT X) 
(ASIN X) 
(ACQS X) 
(ASEC X) 
(ACSC X) 
(ACOT Y X) 

Various trigonometric functions. 



PS:<GLR.FUNCT>EXTFCN.LSP 



PS:<GLR.FUNCT>HYPER.LSP 



(ERROR-FUNCTION X) 

Error function. eps< 1.5E-7. 



(2/sqrt(7r)) Jq^ exp(-t2) dt 



(BESSEL-I N X) 

Modified Bessel function for integer order N. eps < 1.6E-7 

(BESSEL-J N X) 

Bessel function for integer order N. eps < 5E-8. 



PS:<GLR.FUNCT>EXTFCN.LSP 



PS : <GLR . FUNCT>BESSEL . LSP 



PS : <GLR . FUNCT>BESSEL . LSP 



PS : <GLR . FUNCT>GAMMA . LSP 



(GAMMA-FUNCTION A) 

Gamma function. 

r(a)=/o^t^-^e-tdt 

(GAMMA-FUNCTION-INCOMPLETE A X) 

Incomplete Gamma function. <= x <= inf, but breaks if X is large because of roundoff error. 





PS : <GLR . FUNCT>GAMMA. LSP 



r(a,x)--- Vt^"le-^dt 



(BETA-FUNCTION A B) 

Beta function. A and B are FLONUM. 



PS : <GLR . FUNCT>GAMMA . LSP 



^{a, b) = I J t^"^ (1-t)^"^ dt , a>0, b>0 



(BETA-FUNCTION-INCOMPLETE A B X) ps:<glr.funct>gamma.lsp 

Incomplete Beta amction. A and B are FLONUM, but one of A or B must be an integer (eg, 
3.0). 

)g(a,b,x) = /Q'^t^"^(l-t)^'^dt, 0<a,0<b,0< = x< = l 



4 TRIGONOMI-TRIC AND OTHER I-UNCTIONS 



(ELLIPTIC-INTEGRAL-K M) ps:<6lr.funct>ellip.lsp 

Elliptic integral of the first kind, k(M). <= M< 1. cps = 2.0E-8. 

(ELLIPTIC-INTEGRAL-KC M) ps:<glr.funct>ellip.lsp 

Complementary elliptic integral of the first kind, k'(M). k'(M) = k(l-M). 

(ELLIPTIC-INTEGRAL-E M) ps:<glr.funct>ellip.lsp 

Elliptic integral of the second kind, e(M). eps < 2.0E-8. 

/o'^''^sqrt(l-msin2(^))d^ 

(ELLIPTIC-INTEGRAL-EC M) ps:<glr.funct>ellip.lsp 

Complementary elliptic integral of tlie second kind, e'(M). e'(M) = e(l-M). 

(ELLIPTIC-SINE u M) ps:<glr.funct>ellip.lsp 

Elliptic sine function (SN(u,M)) 

(ELLIPTIC-COSINE u M)) ps:<6lr.funct>ellip.lsp 

Elliptic cosine function (CN(u,M)) 



5. Linear Regression 

Say we are trying to fit some (x, y) data to a function that looks like: 

Y'(x) = aO* 1 + al*gl(x) + a2*g2(x) + ... + an*gn(x) 

where tlie a[i] are constants that we are trying to determine, the g[il are linearly independent functions 
tliat we specify, and Y'(x) is the value of the fitted equation. Such a fit is called a linear regression and 
here arc some ftmctions that will do it. 

(FIT FCTN X Y N &optional (MODE 0)) ps:<glr.funct>fit.lsp 

Let X be a one dimensional array of x values (the x values could themselves be vectors). Y is a 
one dimensional array of FLONUM values for the corresponding X input. N is tlie number of terms in 
tlie equation we are fitting (ie, n in tlic equation above). FCl^N is a fiinction of N + 1 arguments that 
evaluates the equation above. That is FCTN looks like: 

(LAMBDA (X AO Al A2 A3 ... AN) 
... ) 

and calculates tJie above equation. FIT returns a list of tJie coefficients Ai as the CDR of its value. 

FIT actually uses Uie procedure RF:GRESS in PS:<GLR.FUNCr>REGRES.LSP, which is 
modeled after <Bevington>. If you want to do something a little more complicated, then ask me about 



SLINI-ARREGRFSSION 



that procedure. 



6. Functional Minimization 

(FMFP FUNCT N X G EST EPS LIMIT) ps:<glr.funct>fmfp.lsp 

Flctchcr-Powcirs functional minimization procedure (adapted from <Kuester> and 
<Luenberger>). This finds a minimum of a multivariate, unconstrained, nonlinear fimction. That is, find 
tlic vector X such tliat F^(X) is a local minimum. The arguments to this function arc: 

N number of independent variables 

X vector[N] of initial variable values 

contains the result vector on exit 
G vector[N] in which to store the gradient 

EST estimate of minimum value of objective function 

EPS test value representing the expected absolute error in movement 

LIMIT maximum number of iterations 
FUNCT user supplied objective function that computes F and the gradient 

(FUNCT N X G) returns F(X) and 

also fills in the vector[N] G with the gradient 

(MARQUARDT N K X Y Z FUNC DERIV B BMIN BMAX BV) 

PS:<GLR.FUNCT>MARQ.LSP 

Marquardt's parameter fitting procedure (adapted from <Kuester>). Given N data points (X[i], 
Y[i]) and a function F with K parameters B[j], find the K parameters B[j] that best fit the data. The fitted 
values are returned in array Z. The arguments are: 

N -- number of data points (FIXNUM) 

K -- number of unknowns (FIXNUM) 

B -- vector[K] of (FLONUM) unknowns 

BMIN -- vector[K] of (FLONUM) minimum values of B[j] 

BMAX -- vector[K] of (FLONUM) maximum values of B[j] 

X -- vector[N] independent variable data points 
(this vector might be a vector of vectors) 

Y -- vector[N] of (FLONUM) dependent variable 

Z -- vector[N] of (FLONUM) computed values of dependent variable 

BV -- vector[K] of (FLONUM) codes 

1.0 -> numerical derivatives 

0.0 -> do not change this unknown 

FUNC -- (FUNC X K B N Z ZOFF) computes the function F (using the 
K parameters in vector B) for each of the X[0] ... X[N-1] 
and puts the respective results into Z[ZOFF+0] ... Z[Z0FF+N-1] 

DERIV-- (not used) 



7 IlYPHRIiOLIC & INVI'RSl- IIYPHRliOl,lC FUNCTIONS 



7. Hyperbolic & Inverse Hyperbolic Functions 

11iis functions use the Common IJsp <Steclc> names, but arc only defined for real arguments. 
They call HXP to compute tlieir results. 

(COSH X) 
(SINH X) 
(TANH X) 
(COTH X) 
(SECH X) 

(CSCH X) PS:<GLR.FUNCT>HYPER.LSP 

Hyperbolic functions. F"LONUM argument and PTONUM result. 

(ACOSH X) 
(ASINH X) 
(ATANH X) 
(ACOTH X) 
(ASECH X) 

(ACSCH X) PS:<GLR.FUNCT>HYPER.LSP 

Inverse hyperbolic functions. FLONUM argument and FLONUM result. 



8. Numerical Integration 

(INTEGRATE-TRAPEZOIDAL F XO XI N) ps:<glr.funct>integr.lsp 

Uses die trapezoidal rule to integrate tlie function F (of one FLONUM argument) from XO to 
XI with N (FIXNUM) iterations. 

(INTEGRATE-SIMPSON F XO XI N) ps:<glr.funct>integr.lsp 

Uses Simpson's rule to integrate the function F (of one FLONUM argument) from XO to XI 
with N (FIXNUM) iterations. 

(INTEGRATE-EULER F XO YO H XI) ps:<glr.funct>runge.lsp 

Uses the forward Bulcr method to integrate F(x, y) from position (XO, YO) to a position (XI, Yl) 
using a step size of H. Yl is the value of IN! HGRATE-EULER. F is a function of X and Y (FLONUM) 
and should return DY/DX for the given position. 

(INTEGRATE-RUNGE-KUTTA F XO YO H XI) ps:<glr.funct>runge.lsp 

This is just like EULER except it uses a fourth order RUNGE-KUri'A metliod instead of 
Eulcr's method. 



8 NUMF.RICAL INTEGRATION 



(INTEGRATE-MULTI-STEP-H F XO YO H XI 

&optiona1 (CD l.OE-8) (CH l.OE-5)) ps:<glr.funct>runge.lsp 

Uses a mulLi-stcp predictor-corrector method to integrate DRRIV (see <Hamming> page 407). 

Integrates from XO to XI using an initial step size of H. When Llie routine estimates that tlie error is 

below CD, the step size is doubled; when the error is above CH, the step size is halved. For 

reasonableness, CD « CH. 



9. Root Finding Functions 

(BISECTION-ROOT-SEARCH F XO XI) 

Uses bisection search to find a zero of F(X) between XO and XL 



PS:<GLR.FUNCT>ROOT.LSP 



(FALSE-POSITION-SEARCH F XO XI EPS) ps:<glr.funct>root.lsp 

Uses false position search to find a zero of F(X) with initial guesses XO and XI. 

(FALSE-POSITION-CONVERGE F XO XI EPS) ps:<glr.funct>root.lsp 

Uses false position convergence to find an X = F'(X) with initial guesses XO and XL 

(CONVERGE F XO EPS) ps:<glr.funct>root.lsp 

Find an X = F(X) with initial guess of XO. Uses Wegstein's method. 



10. Probability and Statistics 

The following functions arc useful in probability and statistics. There are functions for 
computing probability density functions p(x), computing cumulative distributions P(x), and generating 
random numbers from a particular distribution. Most of tlie approximations come from <Abramowitz>. 
Most of the generators use the ratio of uniform deviates method <Kinderman>. 



(UNIFORM-DENSITY X) 
(UNIFORM-CUMULATIVE X) 
(UNIFORM-RANDOM-NUMBER) 

Probability functions for Uie uniform distribution. 

p(x) = l 0<=x<=l 



PS:<GLR.FUNCT>STATIS.LSP 



(NORMAL-DENSITY X) 
(NORMAL-CUMULATIVE X) 

Probability functions for nomial (Gaussian) distribution. 

p(x) = (l/sqrt(2 tt)) exp (-x^/2) 



PS : <GLR . FUNCT>STATIS . LSP 



10 10 PROUAHIl.ITY AND STATISTICS 



(NORMAL-TAIL ALPHA) ps:<glr.funct>statis.lsp 

Given a probability ALPHA, find the X such that ALPHA = 1-P(X), where P(X) is the 
cumulative distribution. 1 his function provides a one-sided tail. 

(NORMAL-RANDOM-NUMBER) ps:<glr.funct>statis.lsp 

Generates a random number from the unit normal distribution. 

(EXPONENTIAL-DENSITY X LAMBDA) 
(EXPONENTIAL-CUMULATIVE X LAMBDA) 

(EXPONENTIAL-RANDOM-NUMBER LAMBDA) ps:<glr.funct>statis.lsp 

Exponential probability functions. 

p(x, X) = Xexp(-Xx), X>0 

(POISSON-DENSITY N TIME) 
(POISSON-CUMULATIVE N TIME) 

(POISSON-RANDOM-NUMBER TIME) ps:<glr.funct>statis.lsp 

Poisson probability functions. Probability of exactly N (FIXNUM) arrivals witliin TIME 
(FLONUM) for a Poisson process with an average arrival rate (X) of L 

p(n, m) = (m" exp(-m))/n!, m = X t 

(CHI-SQUARE-DENSITY X N) 
(CHI-SQUARE-CUMULATIVE X N) 

(CHI-SQUARE-RANDOM-NUMBER N) ps:<glr.funct>statis.lsp 

The CHFSQUARE distribution with N (FIXNUM) degrees of freedom. 

(T-DENSITY X N) 
(T-CUMULATIVE X N) 

(T-RANDOM-NUMBER N) ps:<glr.funct>statis.lsp 

Student's T distribution with N (FIXNUM) degrees of freedom. N must be even in 
T-CUMULATIVE. 

(T-TWO-SIDED X N) ps:<glr.funct>statis.lsp 

For tlic T distribution, the probability tliat T falls within -X to +X. N must be an even 
FIXNUM. 

(F-DENSITY X M N) 
(F-CUMULATIVE X M N) 

(F-RANDOM-NUMBER M N) ps:<glr.funct>statis.lsp 

Snedecor's F distribution witli M and N (FIXNUM) degrees of freedom. 
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11 



(GAMMA-DENSITY X A) 
(GAMMA-CUMULATIVE X A) 

(GAMMA-RANDOM-NUMBER A) ps:<glr.funct>statis.lsp 

Gamma probability functions with parameter A (FLONUM). 1lic random number generator 
must have A > 1.0. 

p(x,a) = (l/^(a))x^"^e"^ a>0 

(BETA-DENSITY X A B) 
(BETA-CUMULATIVE X A B) 

(BETA-RANDOM-NUMBER A B) ps:<glr.funct>statis.lsp 

Beta probability functions with parameters A and B (FLONUM). The random number 
generator is slow, so keep A and B small (say below 10). 

p(x, a, b) = (l/^(a,b)) x^"^ (l-x)^"l, a>0, b>0 

(CAUCHY-DENSITY X) 
(CAUCHY-CUMULATIVE X) 

(CAUCHY-RANDOM-NUMBER) ps:<glr.funct>statis.lsp 

Probability functions for tlie Cauchy distribution. 

p(x) = 1/(77(1 + x2)) 
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